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A Radio Census of Binary Supermassive Black Holes 
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ABSTRACT 

Using archival VLBI data for 3114 radio-luminous active galactic nuclei, we searched 
for binary supermassive black holes using a radio spectral index mapping technique 
which targets spatially resolved, double radio-emitting nuclei. Only one source was 
detected as a double nucleus. This result is compared with a cosmological merger 
rate model and interpreted in terms of (1) implications for post-merger timescales for 
centralisation of the two black holes, (2) implications for the possibility of "stalled" 
systems, and (3) the relationship of radio activity in nuclei to mergers. Our analysis 
suggests that binary pair evolution of supermassive black holes (both of masses > 
10 s M Q ) spends less than 500 Myr in progression from the merging of galactic stellar 
cores to within the purported stalling radius for supermassive black hole pairs. The 
data show no evidence for an excess of stalled binary systems at small separations. We 
see circumstantial evidence that the relative state of radio emission between paired 
supermassive black holes is correlated within orbital separations of 2.5 kpc. 
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1 INTRODUCTION 

Models of the Universe in which hierarchical merging dom- 
inates the growth of galaxies have strong predictions for 
the presence of binary supermassive black holes (SMBHs ) 
at galaxy centres (e.g. IVolonteri. Haardt. fc Madaul [2003). 
As galaxies containing such massive black holes collide, 
the central black holes are expected to inspiral and form a 
bound interacting system which will hav e a significant im- 
pact on the central galactic environment l|Merrittll2006l ). In 
the time preceding their coalescence, these binary sources 
will be among the brightest of gravitational radiation 
sources detectable in the very low frequency ( nHz to /xHz) 
spect rum by pulsar timing experiments (e.g. Ijenet et al.l 
2006). Such experiments are sensitive to both individual 
binary SMBHs with sub-parsec orbits and the stochastic 
gravitational wave background made up of the collective 
signa l from binaries of mass ~10 6 - 1O 9 M (|Detweilerl 
1 19791 ). In the mid-frequency spectrum (10~ 4 Hz to 0.1 Hz), 
the space-based Laser Interferometer Space Antenna 
(LISA] 1 ], is sensitive to gravitational waves from individual 
binary SMBH sytems of mass < 10 7 Mq in the later 
stages of their inspiral, coalescence, and post-coalescence 
ringdown. Predictions of the strength of this astrophys- 
ical gravitational wave background and the expected 
detection rates for LISA and pulsar timing of individual 



* Email: sburke@astro.swin.edu.au 
1 See http://lisa.nasa.gov 



black holes are based on paramet erisations of merger 
rates and the SMBH p opula tio n IIRaiagopal fc Romanil 
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Enoki et al.l 120041: ISesana. Vecchio. fc Colacinol 
Sesana. Vecchio. fc Volonteril 2009), giving these 
tational wave detectors the potential to provide 
insights into the processes of galaxy formation. 

In theoretical treatments, however, the evolution of 
post-merger central supermassive black holes still holds sig- 
nificant uncertainties. The steps in binary b lack hole forma- 
tion a nd evolution were first laid out by iBegelman et al.l 
(1980). After a galaxy pair becomes virially bound, dy- 
namical friction and violent relaxation drive the stellar 
cores containing the massive black holes (of masses mi and 
m.2;rni > 1712) to the centre of the merger remnant on 
roughly a dynamical friction time. When the cores have 
merged, the black holes independently experience dynam- 
ical friction against the merged core's stellar environment. 
This further centralises each black hole, leading the system 
to become a bound pair at a separation where the enclosed 
stellar mass equals the total binary mass. If the stellar core 
is modelled as a single isothermal sphere with stellar density 
p*(r) = a 2 /(2ivGr 2 ), this occurs when the black holes reach 
a orbital semi-major axis of 
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where a v is the velocity dispersion of the merger remnant. 
The length of time spent by black holes to reach this stage 



2 S. Burke-Spolaor 



follows the Chandrasekhar timescale for dynamical friction, 
and is limited by the inspiral time scale for the less massive 
black hole (cf. lLacev fc Coldll993l ): 



tea = 1.654 



Gm,2lnA 



(2) 



where we have assumed circular black hole orbits, and r is 
the orbital radius of the less massive black hole from the 
centre of the galactic potential. In such a system, we can 
assume a value for the Coulomb logarithm InA ~ 5 (a non- 
circular orbit with ellipticity e = 0.5 will add a factor of 
roughly 0.5 t o the timescale and the loga rithm will be 2 < 
InA < 3; e.g. iGualandris fc Merrlttll2008l ). The black holes 
will form a binary when r — at>m- 

The orbital evolution of the black holes proceeds fur- 
ther due to dynamical friction, however as the binary tight- 
ens and increases in velocity, dynamical friction becomes an 
ineffective means of energy transfer and the rate of inspiral 
slows. In this intermediate inspiral stage, there are signifi- 
cant uncertainties in both inspiral mechanism and timescale. 
Relaxing stars in radial orbits enterin g the "loss cone" of the 
massive binary (|Frank fc Reesi [19761 will undergo 3-body 
interactions with the SMBH binary system and be ejected, 
carrying away angular momentum and further shrinking the 
binary orbit. However, without an additional mechanism for 
energy transfer or a means to efficiently refill the loss cone, 
the binary inspiral halts as the loss cone is emptied. Recog- 
nising the long timescal e that loss cone r e -popu lation might 
take for binary systems, iBeeelman et alJ l|l980l ) suggest the 
possibility of gas ejection from the system or gas accretion 
onto the larger black hole, which will cause shrinkage in the 
binary orbit to conserve angu lar momentum (see also, e.g., 
iMerritt fc Milosavlievid [2005). If some intermediate process 
is able to sufficiently shrink the binary orbit, gravitational 
radiation will cause a binary with ellipticity e = to coalesce 
in a timescale 



oc" 



256G 3 mim2(mi + 7712) 



(3) 



N-body simulations and semi-analytical models have 
manifested various inte rmediate inspiral pro cesses (see, e.g., 
the extensive review of lColpi fc Dott1l2009h . The "last par- 
sec problem" is the name given to the hurdle encountered by 
nearly all merger models in which the intermediate inspiral 
mechanisms have timescales sufficiently long that they are 
unable to bring the binary to a regime in which the emission 
of gravitational radiation can drive the black holes to coa- 
lescence in less than a Hubble time. If this is the case, there 
may be many instances of "stalled" binary objects which 
spend a large portion of the ir lifetim es with orbital radii 
within the range 0.01-10 pc (Yu 2002 ). A fu n ction al form 
for the stalling radius is estimated bv lMerrittl (|2006l ) : 
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If systems are unable to find the fuel to reach coalescence 
in less than the age of the Universe, then they may appear 
to stall indefinitely at orbits greater (thus GW frequencies 
lower, and amplitudes generally smaller) than those in the 
range detectable by either pulsar timing arrays or LISA. 
This would have dramatic consequences for GW detectors 
and for hierarchical formation models alike. 



Since few small-orbit binary systems are known and 
electromagnetic tracers of post-merger galaxy cores are hard 
to identify, it has been difficult to study the post-merger 
dynamics of binary black hole systems. With only two con- 
firmed sub-kpc double b lack holes serendipito usly identified, 
fNG C 6240 at ~1 kpc, Ikomossa et all 120031 ; 0402+379 at 
7 pc. lRodriguez et al'll2006l ) there is little possibility to assess 
the physical mechanism responsible for driving binary sys- 
tems into the phase when gravitational radiation dominates 
the binary inspiral. Observationally, there are several ways 
to identify galaxies wi th binary nucle i, which are reviewed 
briefly below (see also iKomoss al l2006l for a good overview ) . 

Pairs of galaxies/quasars with small angular separa- 
tions and similar redshifts, as well as galaxies which show 
disturbed structure and tidal tails in the optical and in- 
frared, can indicate interaction betwee n widely separated 
stellar cores containing SMBHs (e.g. I Myers et al.l 120081 : 
iDe Propris et al.ll2007l ; IWen et al.ll2008D . Secondary signals 
that have been modelled as a consequence of a close or bound 
binary SMB H include quasi-perio dic oscillations in flux den- 
sity fOJ287. IValtonen et al.lll988l ), abn ormal jet morphology 
such as helices or X-shaped radio jets |Saripalli et al.l [20081 : 
lLal et all I in press! ) , spatial oscillations of an active galac- 
tic nucleus (3C66B, ISudou et al] l2003h . and sources with 
double or offset broad/narrow line regions, implying active 
galactic nuclei (AGN) moving at high relative velocities to 
each other or to the host galaxy (e.g. Peterson et al. 1987 : 
Comerford et al.ll200Sl ; lBoroson fc Laueij|2009l : lDecarli et all 
20101 ). In most of these cases, other explanations have 



better suited the observed characteristics of each source. 
However, spatially resolved systems such as 3C75, which 
shows two jet-emittin g radio AGN cores, and NGC 6240 
l|Komossa et alJ 20031) and o ther double X-ray emitting 
AGN fe.g. lWang fc Gaolboiol ). are straightforward to iden- 
tify and can provide a direct detection of a double SMBH 
system. 

The core regions of radio-luminous AGN have been ob- 
served to show a distinct radio continuum spectrum which 
typically peaks at radio frequencies above ~lGHz, giving 
core components a characteristically flat two-point spectral 
index (a > — 0.5, where radio flux F v oc v a ) at GHz fre- 
quencies (e.g. ISlee et al.|[l99i ). This feature distinguishes 
the region most directly associated with the host black hole 
from the outer features such as larger scale jets or hotspots, 
which show power-law spectra generally at or steeper than 
a = —0.7. Taking advantage of these unique spectral signa- 
tures of the core regions of radio-emitting AGN at high radio 
frequency (1 to 20 GHz) and the sub-milliarcsecond precision 
of observations using radio Very Long Baseline Interferom- 
etry (VLBI), it is possible to spatially resolve and identify 
AGN containing multiple cores using multifrequency giga- 
hertz spectral imaging. This property has also been utilized 
to locate gravitationally lensed compact radio sources; if an 
AGN is lensed, it will also appear as a spatially resolved pair 
with equivalent radio continuum spectra (note that there is 
dire ct overlap in the targe t of such such searches; because 
e. g. IWilkinson et al.ll200ll found no lenses in 300 VLBI im- 
ages, this implies no SMBH binaries were found, either). 

In the case of a genuine physical pair, resolutions of 
1 milliarcsecond are able to discern components with a pro- 
jected spacing of above ~8.5pc at all redshifts, and down 
to sub-pc resolutions for the nearest galaxies. This makes 
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Figure 1. The redshift distribution of 1575 objects in our sample. 
The closest source is at z = 0.000113, while the highest redshift 
source is a quasar at z = 4.715. 



VLBI observations well-suited to explore double supermas- 
sive black hole systems even at very small projected sepa- 
rations, and within post-merger galaxies that are otherwise 
unidentifiable as double or disturbed systems. Identification 
of double black holes in post-merger systems and binary 
black holes in the intermediate stage of inspiral will aid in 
showing evidence for stalling and may provide statistical es- 
timates of the rate of SMBH inspiral at various phases. 

This paper reports on the first systematic search of 
a large number of radio sources for spatially resolved bi- 
nary radio AGN using archival VLBI data. We present a 
brief description of the archival data sample in Sj2] and sum- 
marise our search technique in Sj3] The candidates resulting 
from this search are assessed in f|4] We outline an interpre- 
tive framework for making statistical estimates of inspiral 
timescales based on our detections in ij5l and give the results 
of this analysis in ij6] In Sj7]we discuss the analysis in terms 
of implications for black hole inspiral rates, stalling, radio 
core emission in the merging process, and implications for 
the pulsar-timing-detectable gravitational wave signal from 
SMBH binaries. 



2 THE ARCHIVAL DATA SAMPLE 

The data used in this search were extracted from the VLBI 
archives maintained online as the Goddard Space Flight 
Center astrometric and geodetic catalogues^ The catalogue 
version used was "VLBI global solution 2008a_astro," which 
contains a large volume of calibrated VLBI data from 1980 
to ea rly 2008, including data from the VLBA calibrator sur- 
veys jBeasley et al.ll2002l;lFomalont et al.ll2003l;IPetroy et all 



120051 . 120061 ; iKovalev et all 120071 ; iPetrov et all l20Qgft . The 
catalogue contained observations unevenly sampled across 
2.2 GHz to 43 GHz for 4169 radio sources. For our search 
technique, we are limited to sources observed at two or more 
frequencies, and to avoid potentially problematic intrinsic 

2 http://lacerta.gsfc.nasa.gov/vlbi/solutions/ — Please see web- 
site and references therein for a complete description of original 
data sets 
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15-24 190 
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Table 1. Frequency band pairs usable in our search; all frequen- 
cies are in units of gigahertz. The frequencies above refer to the 
rough observing band only. In total there were 3982 frequency 
pairs for 3114 sources. 



source variability (see 331, we rejected frequency pairs which 
were observed at times separated by more than 35 days. Our 
sample was further culled because of software incompatibil- 
ity for some data, and for the remaining dataset, observa- 
tions were rejected if the images revealed data which con- 
tained no clear source above the noise. Three sources were 
removed because a literature searched revealed them to be 
Galactic HII regions or stars. In total, we were able to search 
data for 3114 sources. 

Approximately half of the sources in our sample had 
redshift information available. The distribution for the 1575 
sources of known redshift is shown in Fig. [1] the source 
types and redshifts were gathered from the NASA Extra- 
galactic Database^ The mean redshifts of our sample are 
z g — 0.402, z q = 1.401, and z to t = 1.226 for galaxies, 
quasars, and all the sources with known redshift, respec- 
tively. 



3 VLBI SEARCH FOR DOUBLE AGN 

Our technique for identifying binary active galactic nuclei 
using VLBI is essentially a search for objects which show 
more than one flat-spectrum component. It requires mor- 
phological analysis and spectral index imaging across two 
or more frequencies to distinguish the nucleus of an object 
from its other components, such as complex extended struc- 
ture or bright, unresolved hotspots along a jet. These other 
source components, as was previously noted, show system- 
atically declining synchrotron spectra at a < —0.7. The 
technique and the use of archival VLBI data require sev- 
eral precautions, primarily because both the accuracy of 
frequency-dependent flux density measurements and mor- 
phological fidelity of the radio image are important to this 
method. 

We therefore take special precaution in considering 
AGN variability and coverage of the spatial-frequency (u—v) 
plane for the data sets. The variability exhibited by some ra- 
dio AGN typically comes in two different flavours. Both par- 
ticularly affect compact regions of emission with high bright- 
ness temperature, such as radio cores. The first is intra-day 
variability (IDV) due to propagation effects in the inter- 
stellar medium. This signal flickering can occur on sub-day 
timescales and exhibits a significant modulation index at 
frequencies below 5 G Hz, causing fl u ctuations on avera ge at 
the 1-10% level (e.g. I Walker! 1 19981 ; lLovell et alJl200Sh . In- 
trinsic variability in radio AGN cores is observed most sig- 
nificantly at frequencies above 5 GHz, and occurs on a wide 
range of timescales. Some radio AGN exhibit slow variations 
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RA offset (arcsec; J3000) RA offset (arcsec; J300D) RA offset (arosec; J3000) 

(a) J0356+6043 (b) J1829+4844 (c) J1823+7938 



Figure 2. Two-frequency spectral index maps for sources tagged as (a) pointlike, (b) resolved, and (c) multiple fiat-spectrum component 
sources. All maps show matched-resolution data. The color scale represents the spectral index calculated from two input maps of differing 
frequency. The color scale is linear in greyscale from —1 < a < —0.5, and linear in color from —0.5 < a < 1. The gray contours show the 
shape of the source at the lower frequency used to calculate the spectral index map, while blue contours trace the source morphology at 
the higher observing frequency used to calculate the spectral index map. These sources were randomly chosen from the source list for 
each category; all show matched-resolution data between 2.3 GHz (gray contours) and 8.4 GHz (blue contours). Contours in these plots 
are set at 2, 25, and 50 % of the peak flux at each frequency. 



on timescales of years. The existence of such variations place 
limitations on the allowable temporal spacing between ob- 
servations at different frequencies. Instantaneous multifre- 
quency spectra will not be influenced by either effect, while 
observations spaced at less than the characteristic timescales 
of long-term AGN variability will allay the effects of a wan- 
dering high-frequency signal. For this reason we used only 
frequency pairs with a time separation of less than 35 days. 
Out of the 3982 usable frequency pairs in our data, 3447 
were taken simultaneously, leaving 535 with an average dif- 
ference in observation time of about 11 days. 

All data processing was done using the Miriad soft- 
ware package (|Sault. Teuben. fc Wrightlll993h . For each fre- 
quency pair, we inverted the u — v data, tapering the maxima 
of u and v to give a circular synthesised beam of the same 
size at the two frequencies. A manual inspection of u — v cov- 
erage for the pairs ensured that the sources were sampled 
enough to provide an accurate representation of the source 
morphology at the tapered resolution. Images were cleaned, 
restored with the synthesised Gaussian beam of the tapered 
resolution to make images of 512x512 pixels, allowing a typi- 
cal field of view 1.2 arcseconds to a side (though varying with 
the resolution of the observation) . An estimate of the image 
noise level at each frequency was determined using an off- 
source region of each image. The number of frequency pairs 
we were able to use is detailed in Table Q] 

The data were imaged in three ways for manual candi- 
date searching. First, for each resolution- matched frequency 
pair, the full-field-of-view images at low and high frequency 
were plotted in greyscale with logarithmic contours begin- 
ning at three times the noise level of each image. Enhanced 
plots of the central image regions were made in a similar 
fashion, and a spectral index map of the inner regions was 
computed, calculating the spectral index pixel-by-pixel at 
any point in the image where the flux at both frequency 
bands was greater than three times the image noise. Spec- 
tral index images were superimposed with contours for each 
respective frequency, and the spectral index value was plot- 
ted as a colour image with a sharp desaturation break at 



a = —0.5 to enhance the appearance of flat-spectrum source 
components. 

We then manually assessed each source, dividing 
sources into pointlike, resolved/extended, and multiple flat- 
spectrum component sources. Representive sources of each 
category are shown in Figure [2] Sources with multiple flat- 
spectrum components were considered preliminary candi- 
dates. In cases where the sources in multi-frequency images 
appeared displaced, a manual best-fit alignment was done 
and the source was reassessed. We considered a component 
to be "flat spectrum" if its two-point spectral index was 
a > —0.6 at any frequency pair, or if structure was visible 
at high frequency which was not detected at low frequency. 
As previously noted, sources which contained no detection 
above a 4<r noise level were removed from the sample. 

The nature of preliminary candidates was then assessed 
by imaging all available observations of each candidate. 
Sources were then rejected from the candidate list if they 
satisfied any of the following: 

- In the case of discrete compact fiat spectrum compo- 
nents, if suspected components showed spectral steepening 
at later observation epochs (implying a false detection due 
to a young, evolving jet component). 

- If an analysis of the time-dependent movement of sus- 
pected flat-spectrum source components revealed the com- 
ponents) as members of a steady outward (in some cases, 
superluminal) flow along a jet axis. 

- If the full resolution image revealed beam sidelobes as 
the source of a spurious, flat-spectrum component. 

Our search returned 12 candidates containing genuine 
multiple flat-spectrum components as discussed in more de- 
tail below. 



4 CANDIDATES AND SEARCH RESULTS 

Of the 3114 sources in our search, 68 percent were found 
to be pointlike at our frequency-matched resolution, while 
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(a) Candidate J0541-0541 (b) Candidate J1658+0741 (c) Candidate J2246-1206 



Figure 3. Example two-frequency spectral index and contour maps for three of the initial candidates described in £14,11 Coloring of the 
mapped data and contours are as in Fig. [2] Contours are set at 4, 16, 64, and 256 times the RMS image noise drawn from the respective 
frequency's input image. 



the remaining sources were identified as resolved or other- 
wise complex systems with only one flat spectrum compo- 
nent. The number of sources which exhibited multiple flat- 
spectrum components comprised less than 4 percent of the 
sample, and after the assessment of preliminary candidates, 
there remained only 12 candidate binary sources. As an ex- 
ample of our data quality for candidates, spectral and two- 
frequency contour images for three of the candidates be- 
low are shown in Figure One c andidate was the binar y 
black hole candidate published bv lRodriguez et"aL | (|200fJ ): 
0402+379; an investigation of this source is available in their 
publication. We found that the remainder of our candidates 
appear to contain regions with active jet injection, which 
effectively gives rise to a population of contaminants for the 
radio spectral index mapping technique, in some cases re- 
quiring further data and analysis to distinguish such a source 
from a genuine binary black hole. Below we use our data 
and that available from the literature to further investigate 
each source and gather evidence for the presence of such 
processes, and not a binary black hole, in the sources. 

Ultimately, we report the detection of no new binary 
systems. We have searched a sample of considerable size, 
finding only one spatially resolved binary black hole out of 
3114 radio-bright AGN. Therefore, within the limitations of 
our search and data, we find 0.032% of sources in our sample 
in which a binary radio AGN was detected. The implications 
of this are discussed in later sections. 



4.1 Candidates with active relativistic injection 

Eleven of twelve candidates were rejected after the exami- 
nation of all available data (in our data, the literature, or 
available onlin^B) revealed that the multiple fiat-spectrum 
components were part of a jet. These sources could be dis- 
criminated from a binary SMBH by proper (superluminal) 
motion, or by the resolving of the emission in higher reso- 
lution maps and alignment of the components with a more 
extensive outflow: 



We made use of the MOJAVE database at 



JO 11 1+3906: The dynamical age o f the source (370 yr, 
lOwsianik. Conway, fc Polatidis|[l99Sl . and the extendedness 
and alignment of the eastern component along the larger 
source structure indicate that it is a young jet. 
J0541-0541: The southwestern flat spectrum component 
contains no compact emission at the high-resolution 8.4 GHz 
observation. This source is shown in Figure [3(a) | 
J0741+3112: T he persistent total intensity and polarisation 
observations bv lLister et al.l l|2009h clearly demonstrate that 
this source has a continuous outflow which is fuelling the flat 
spectrum hot spot appare n t in o ur maps. 
Jl 147+3501 : [G iovanninll (|l999h observe superluminal mo- 
tion of the w estern flat-spectrum components . 
J1223+8040: \Po\lack, Taylor, fc Zavalal (|2003t ) indicate that 
there is a jet traversing our flat-spectrum components. 
J 1347 +1211: The multiple flat-spectrum components are all 
members of an active superluminal outflo w, as shown by the 
five year modelling of iLister et al.l l|2009T ). 
J1459+7140: The southern component is fully resolved at 
our highest resolutions, and aligned with a small-scale jet 
emerging from the northern flat-spectrum component. 
J1658+0741: The weaker flat-spectrum component of this 
source is resolved in our non-tapered images and is aligned 



with the northern jet, visible in Figure 3(b) 
J 1823+7938: The resolved structure of the easternmost 
component in our d ata and the proper motion analysis of 
iBritzen et alj (|200&t ) indicate that it is not a binary radio 
nucleus. This candidate is shown in Fig. 2(c). 
J2246-1206: We conclude from the motion analysis of 
([Lister et al l [2009) that the northern component is a jet 



member. This candidate is displa yed in Figure 3(c) 
J2253+1608: ILister et al] |2009r i show the western flat- 
spectrum component to be a jet member. 



http://www.physics.purdue.edu/MOJAVE which is maintained 
by the MOJAVE team ILister et all <2009h . 



5 FRAMEWORK FOR INTERPRETATION 

The use of radio data as a black hole indicator puts immedi- 
ate constraints on both the physical properties and number 
of sources we will be sensitive to. We adopt the following 
expression for the number of sources expected to be found 
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in our search: 



A^cxp = ^2 fbbh(zi,mii m ) ■ Pi 



(5) 



Here, N s is the number of objects that have been searched. 
It is assumed that a target's redshift z is known (true for 
1575 of our objects), and that it contains at least one de- 
tectable black hole of mass M, > mii m P represents the 
probability that we would detect a second black hole if the 
system is a binary, being the probability that the secondary 
black hole is both radio-emitting, and bright enough for a 
successful detection in our images. This term is explored in 
Sec. 15.11 below. The fraction of galaxies, /bbh, containing a 
SMBH binary system at a redshift z' to which we are sen- 
sitive will be the cumulative number of galaxy mergers over 
the epoch at z' until the epoch at z" [t z — t v i s ] , where t z is 
the age of the universe at redshift z' , and t v i s is the dura- 
tion that a signature of a double black hole in the merging 
system is detectable in our search. Thus, if a binary stalls at 
a separation detectable by our technique (bounded by our 
observing resolution and the field of view of our images; see 
£|5.3[) . galaxies that began to merge when the Universe was 
the age of z" will still be occupied by a resolvable binary 
at z . The occupation fraction for binaries at a resolvable 
separation is then given by: 



/bbh (2, TOlim) = 



/*/ A r mrg(2,mi im ) d,Z 

iV ga i(2',mii m ) 



(6) 



where -/V mrg is the number of virialized galaxy pairs with 
SMBHs of masses above mnm at the given redshift, while 
iVgai is the total number of galaxies at z with SMBH masses 
greater than mii m . The timescale of inspiral over the scales 
to which our data are senstive (t v is) can thus be estimated 
with a prediction or measurement of the redshift- dependent 
merger rate of galaxies containing supermassive black holes. 

Whether we can identify a binary depends on factors 
explored in the sections below. The timescale over which we 
can resolve a binary is determined by our spatial sensitivity, 
the black hole masses and binary mass ratio; the relative 
time spent at various separations in a binary's inspiral de- 
pend on black hole mass. If the radio luminosity of the black 
holes somehow depends on the stage of the merger (for in- 
stance if the radio mechanism shuts down as the galaxies 
begin interacting, or ignites only after a binary coalescence), 
an additional factor needs to be considered in Eq. [5] How- 
ever, because the evidence for a relationship between radio 
fuelling and merger events has until now been tenuous and 
estimates for an AGN timescale are so far inconclusive, we 
limit this uncertainty by employing two scenarios below that 
represent the most optimistic and the most pessimistic ex- 
pectations for the state of radio luminosity in the two black 
holes. 



5 Throughout our analysis we use a value mi; m = 10 s Mq; here 
we make use of the observed property that radio-bright quasars 
have a distinct mass distribution: (log( M./M a )) = 8.89 ± 0.0 2 
with a lower cutoff around 10 8 Mq (e. g. iMcLure &: Jarvi 12004) . 
We probe only this most massive end of the black hole population, 
and for a double detection to be made, the second black hole must 
exceed a limiting mass mi; m ^ 10 8 Mq. 



5.1 Supermassive Black Holes and Radio 
Emission 

All our targets are radio sources detectable on parsec scales, 
and it follows that all contain at least one radio-emitting 
supermassive black hole. If a target is a binary system with 
a resolvable separation between the two black holes, due 
to our inspection techniques a second radio-emitting SMBH 
will have only been flagged as a candidate if its flux den- 
sity exceeds four times the root-mean-squared noise (o" rm s) 
in its image. This 4 cr rms flux limit corresponds to a luminos- 
ity sensitivity limit of Lu m (0> ms , z). To predict the expected 
number of binary AGN to be found in our search and es- 
timate a timescale for binary SMBH inspiral, we must first 
quantify the likelihood that the radio emission of a second 
black hole will exceed this luminosity limit. 

The main unknown contributing to the discussion is 
how strongly the galactic and intergalactic environment re- 
late to the production of radio AGN activity around the cen- 
tral black hole: that is, do the merger and mixing of gaseous, 
dusty, and stellar environments contribute strongly the pro- 
duction of radio activity? Below we assign two scenarios to 
bound assumptions at either end of extremes: 1) Pessimistic; 
the probability that a SMBH paired with a radio-emitting 
AGN is itself radio-emitting is no greater than the proba- 
bility that a non-binary SMBH at the center of a massive 
galaxy would be radio-emitting, and 2) Optimistic; given a 
black hole that exists in a shared environment with a radio- 
emitting AGN, the probability of radio emission in the sec- 
ond black hole is unity. 

Given that the black holes in the progenitors may not 
encounter significantly mixed environments until the binary 
has a relatively small separation, it is also possible that ra- 
dio ignition is radially dependent. Regardless, we stress that 
for the third stage of inspiral in which theoretically the bi- 
nary may stall, the stellar cores around the black holes have 
merged and the black holes will have shared many orbits 
within one another's sphere of influence and in a common 
environment. 



5.1.1 Pessimistic scenario 

We can make an estimate of the fewest number of bi- 
nary sources we would expect to see by setting P equal 
to the integrated bivariate luminosity function of galax- 
ies, ip{M, >Ln m ), where M represents the absolute opti- 
cal/infrared magnitude of a galaxy. In this case we assume 
that mergers have no influence on the radio-active state 
of a black hole; that is, that the probability that the sec- 
ond black hole in a binary would have a detectable lumi- 
nosity L > Ly lm is equivalent to the radio active fraction 
for non-merging galaxies selected from a random pool with 
M, ^ miim. We calculate P for each source using a log-linear 
in terpolation of the integr ated bivariate luminosity function 
of iMauch fc Sadlerl (|2007l ) for radio AGNs in galaxies with 
—24 > Mk ^ —25, taking such galaxies to be represen- 
tative of our tar g ets' h ost galaxies. As roughly 10% of the 
IMauch fc Sadler] (|2007|) sample have flat spectra (based on 
the measurements of Mason et al.l |2009h ). we expect that 
approximately 10% of their sources will be compact down 
to VLBI scales. To account for this we divide the space den- 
sity given in their radio luminosity function by a factor of 
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10. For our 1575 sources with a known redshift, we find that 
in this scenario P has a mean and standard deviation of 
0.4% and 1%, respectively. 



5.1.2 Optimistic scenario 

In the optimistic scenario, we consider that a SMBH in the 
presence of a radio AGN will itself be a radio- loud AGN. 
This scenario implies that when the two black holes share a 
common galactic environment, both black holes will be radio 
luminous if the conditions are fit for radio ignition; that is, 
both objects will be either radio-loud or not. We take as a 
given that one black hole is a detectable radio AGN in each 
of our images. To determine the detectability of a second 
SMBH in such a system (the probability that its observed 
luminosity is greater than Lii m ), we would like to avoid the 
explicit use of e. g. relativistic beaming models, correlations 
between black hole mass and radio luminosity, and the levels 
of flux resolved by our baselines — the parameters of which 
are uncertain — by using an estimate derived from observa- 
tions in which such effects are inherent. 

We therefore take the secondary SMBHs in such sys- 
tems to have a radio luminosity distribution equal to the 
luminosity distribution of solitary, flat-spectrum radio black 
holes, which will implicitly include beaming and mass- 
correlation effects. The use of the flat-spectrum popula- 
tion assumes that the emission observed in flat-spectrum 
radio sources is largely unresolved, and therefore scal- 
able to the resolutions reached by our VLBI data. We 
adopt a probability density function for secondary black 
hole radio luminosities, p(z, L) dL, based on the redshift- 
de pendent, flat-spectr u m "pu re luminosity evolution model" 
of iDunlop fc Peacock! (|l99dl ). We may then determine the 
probability that for a given object in the sample, a second 
black hole will have a luminosity above our detection thresh- 
old, Lum- The probability that we would have detected the 
second black hole of a binary AGN is therefore: 



P=— p(z,L)dL. 

Pn J L lim 



(7) 



— p(z, L) dL, sets the 
of radio-loud AGN. 



The normalisation factor, p 

range of radio luminosities of radio-loud AGN. We set 
the lower luminosity bound of radio A GN at Lrl = 
10 23 W Hz" 1 sr -1 , equal to that given by IPadovani et all 
(2009). Although the exact value of Lrl debatable, the av- 
erage P value is not highly sensitive to the value of p n if 
Lrl is decreased. For images in which the Ln m < Lrl, the 
probability of detection is set to one. For the 1575 known- 
redshift sources in this scenario, P has a mean and standard 
deviation of 57.9% and 23.7%, respectively. 



5.2 Host Properties and Limiting Radii 

The predicted binary formation radii and stalling radii mark 
key points in possible inspiral rate changes for a post-merger 
system. In order to determine what stage of inspiral we have 
probed with our spatial sensitivity limits, we require an esti- 
mate of the predicted typical abin and OstaU for our sample. 
As mergers between less massive galaxies are expected to 
be more common than high mass mergers, for this study we 
expect the least massive (and most common) merger probed 
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Figure 4. The upper and lower bounds of the average linear 
scales searched as a function of redshift. The error bars represent 
the standard deviation of (D m i n ) and (D max ). 



100 



10 



-9= 1 



0.1 



0.01 



0.001 











1 >) 




































Slll^x 




















A *&> 


























a stal 
a biq 




- . I 









0.0001 



0.001 



0.01 0.1 
redshift 



10 



Figure 5. The minimum resolvable linear scale for each of the 
sources of known redshift in our sample. Overplotted are the pre- 
dicted black hole binary formation radius and stalling radius for 
two 10 s Mq black holes. Note that for a binary in a circular or- 
bit, these radii correspond to a maximum projected separation of 
twice each value. 



in our sample to involve two black holes of 10 Mq. Relevant 
progenitor galaxy properties can be calculated to estimate 
the characteristic scales (e.g. stalling radius, binary forma- 
tion radius) that our study probes. Using two black holes 
of mass 10 s Mq, each progenitor's velocity dispersion, a v ~ 
187kms _1 , and bulge mass, M hoat ~ 6.63 x 10 10 Mp), can 
be est imated by the M. — a v relation of Fcrrarese fc Merrifrtj 
(2000) and the M. — Mbulge relation of lHaring fc Rixl l|2004h . 



respectively. If the black holes do not accrete a significant 
fraction of their mass during inspiral, the post-merger galaxy 
will contain a black hole of 2 x 10 s Mq and an implied ve- 
locity dispersion of a v ~ 216kms _1 . The binary formation 
radius for our typical system is thus abm = 27.7 pc, while 
the stalling radius is a sta n = 2.8 pc (equations [T] [4| . In the 
earliest stage of merger, the virially bound galaxy pair will 
have a projected separation < 7kpc. We note that while 
these radii are the characteristic values for the sources in 
the sample, more massive (and less common) mergers be- 
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tween 10 Mq black holes will have larger radii by a factor 
of ~10. 



5.3 Spatial Limits of the Search 

Our technique resolves a large range of spatial scales, over 
which an inspiral rate will undergo several changes. To es- 
timate post-merger black hole evolutionary timescales we 
must determine the projected spatial separation range that 
our data set is sensitive to. For the 1575 sources of known 
redshift, the average maximum and minimum projected lin- 
ear spatial sensitivities ((D mlx ), (D m i u ), both in units of par- 
sec) as a function of redshift are plotted in Figure fj] The 
sources are binned in increments matching the redshift steps 
of the Millennium Simulation snapshots (see §5.4\ each step 
spans an average of 300 Myr) . 

The observed scales are derived as follows: the min- 
imum and maximum linear projected separation between 
two objects that we are sensitive to are set by our angu- 
lar resolution, 9 r , and half the field of view of our search 
maps, given that one AGN will lie at the pointing cen- 
tre of the observation. Using the redshift of each source 
and a flat-universe cosmology with Ho = 72kms _1 Mpc~ , 
fim = 0.27, and Q,a = 0.73, we calculate d z /8 z , the distance 
(in pc) per milliarcsecond at that redshift on the plane of 
the sky. Our smallest resolved size for that source is then 
D m i n = dz6 r /6z parsecs (D m i n for individual sources are 
plotted in Figure [5]). Because we searched images of 512 
pixels to a side and the synthesised restoring beam is 3 pix- 
els in breadth, the maximum linear scale we are sensitive to 



512 8 r d z 
fmax = — - — • -p- parsecs . 



(8) 



This value represents the largest scale size to which the 
search of each source is complete. For sources with known 
redshift, the sample averages were (9 r ) = 6.5 mas, {D m in) = 
40 ± 25 pc and (D max ) = 3415 ± 2133 pc. 

The limited spatial separation window over which we 
have searched each source for a double AGN and the poten- 
tial for drastic changes in inspiral rate at various stages of 
the binary evolution compel us to divide our measurement 
of inspiral rates into three stages. We do so by segregat- 
ing the sources by their spatial resolution into groups which 
represent sensitivity to three stages of merger evolution: the 
stage at which the stellar cores have not yet merged (merger 
stage 1), the stage leading to binary formation in which the 
cores merge and the binary forms (stage 2), and the final 
stage to which we can probe spatially, in which the black 
holes are a binary and evolve to the predicted stalling ra- 
dius (stage 3). In each of these groups, we consider that a 
linear resolution of D m i n = ai will allow the detection of 
sources with an orbital semi-major axis of ai roughly 80% 
of the time, assuming systems with a random distribution 
of inclination angles. By setting the constraint on a group 
that -D m in < aio, we ensure that the searched nuclei in that 
group will be sensitive to a binaries with axis a ^ ai for 
> 80% of the pair's orbit. Each group is statistically sensi- 
tive in this way to pairs with an orbital axis of up to (D ma x} 
of the group. 

The three groups are divided at (1) 
27.7 < D m in < 120 parsecs, (2) 2.8 < D min < 27.7 parsecs 
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Figure 6. The redshift-dependent fraction of galaxies containing 
binary or double SMBH black hole systems for various inspiral 
times , based on the merger rate predictions of lDe Lucia fc BlaizoiJ 
for mergers with black holes of masses mi, mi 10 8 Mq. 



(sensitivity down to our limiting awn), and 
(S) Dmin < 2.8 parsecs (giving sensitivity down to the 
minimum a sta ii). Because sources with very high resolution 
observations will not have a large value of D max and thus 
not probe large separations, the source group (3) is not 
included in (2), and in turn these two groups are not 
included in (1). The divisions give a sensitivity to binaries 
at separations over three ranges of scales: 



(1) 120 < a < 4757 parsecs 

(2) 27.7 < a < 2519 parsecs 

(3) 2.8 < a < 450 parsecs 



There were 1035, 497, and 43 sources in groups one, two, 
and three, respectively. 



5.4 Galaxy Merger Rates 

The merger rate of galaxies containing SMBHs given a 
ACDM universe in which merging is the primary mech- 
anism for galaxy growth can b e predicted using r esults 
from the Millennium Simulation (|Springel et al.ll2005T ). The 
Millennium simulation is a large-scale N-body simulation 
which tracked the evolution of dark matter halos in a 
co-moving cubic volu me SOOMpc/i -1 to a sid e (h — 
Hn/lOOkms" 1 Mpc" 1- ). iDe Lucia fc Blaizotl (120071 ) applied 
a semi-analytical prescription to the Millennium simulation 
to track the evolution and merger histories of individual 
galaxies within the Millennium volume. Their catalogue can 
be queried through an online interface^ 

Using the Millennium catalogue, we constructed the 
redshift-dependent distribution of galaxies and mergers in 
which each galaxy contains a supermassive black hole of 
M. > 10 8 Mq . For comparison with the black hole binary 
occupation fraction, we show these values as a pair frac- 
tion, where / mr g = Nmxg/Nggi; here, iVmrg and iV ga i are the 
number of mergers and the total number of galaxies, respec- 
tively, with M m > 10 8 Mq at the given simulation snapshot. 



See http: / /www.g- vo.org/Millennium 
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A "merged" pair in the Millennium Simulation and the De 
Lucia & Blaizot framework is a pair which has become grav- 
itationally bound in the ~300Myr between two adjacent 
snapshots. 

Figure [6] shows the redshift-dependent pair fraction for 
mi im > 1O 8 M and the corresponding black hole binary oc- 
cupation fraction (computed from eq. [5]) for a range of in- 
spiral times. The merger distribution peaks at a redshift of 
about z = 0.7. 



6 LIMITS ON INSPIRAL TIMESCALES 

Following Eq.O for the merger stage groups with no SMBH 
binary detections, we place an upper limit on tvis at 7V cxp = 
0.5. The one source in which a binary black hole was de- 
tected is part of the second merger stage group; for this 
group/region only are we able to make a measurement of 
(rather than put a limit on) tvis, with errors bounding where 
tvis gives a value of 0.5 < iV oxp < 1.5. 

The timescale measurements for the optimistic scenario 
are summarised in Figure which illustrates the number of 
black holes expected to be found in each group as a function 
of tvis, and Table [21 which gives the numerical values of our 
limits. Under the pessimistic set of assumptions, we would 
not expect to see any sources, given that even with t v - m — 
l/Ho, we find 7V cxp = 0.025,0.20, and 0.16 for groups 1, 2, 
and 3, respectively. The implications of this will be discussed 
in the next section. 

In Table[5]we report the timescales measured from both 
our 1575 sources of known redshift, and from the full sample 
with an assumed redshift distribution. While we do not use 
the less reliable "full sample" timescale values in the analysis 
below, we stress the limiting power of the sample if we were 
to have a full sample with measured redshifts, by assuming 
the remaining sources have a redshift distribution equivalent 
to the other 1575. This has the effect of roughly doubling 
iV cxp for each group, and as such works only to lay more 
stringent limits on i v is- 

Finally, although the number of sources in group one 
was the largest, these sources also had a tendency to be 
of higher redshift, with a net effect of fewer expected de- 
tectable binary black holes due to the lower merging fraction 
of SMBH hosts at redshifts z > 1. For stalling studies, the 
observation of low-redshift sources is therefore doubly ben- 
eficial: not only does it allow smaller scales to be resolved 
in the sources due to their proximity, but it affords a higher 
predicted success rate due to the increased merger rate of 
SMBH-hosting galaxies at low redshift. 



7 DISCUSSION 

The question we aim to answer with this search is whether 
there is observational evidence for stalled supermassive bi- 
nary black hole systems: that is, what is the evolution 
timescale down to the stalling radius for a post-merger su- 
permassive black hole pair? For our methods, this question 
is intertwined with how closely correlated the radio emis- 
sion is between two black holes in a pair. Our sample al- 
lows us to address this correlation and inspiral timescales 
over three stages of post-merger SMBH pair evolution. This 
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Figure 7. The cumulative expected number of sources and lim- 
its on the inspiral timescale for each spatial-sensitivity limited 
group using the optimistic scenario f t| 5 . 1 . 2 D in which given one 
radio-active black hole in a binary, the second will also be radio- 
emitting. For groups 1 and 2, the upper limit set on t v - ls is marked 
with an arrow, while for group two the measured value of t v - ls is 
marked by a filled circle, with open circles indicating the error 
range. 

Known redshift Full sample 

Group 1 <0.119 <0.060 

Group 2 0.2861^3 0.143 

Group 3 <0.846 <0.416 

Net <1.25Gyr <619Myr 

Table 2. A summary of the t v - ls measurements (in Gyr) placed 
by our known redshift sources (1575) and the full sample with 
an assumed redshift distribution (3114 total). The values given 
are for the set of optimistic radio luminosity assumptions f M5- Xp . 
Groups 1-3 are as detailed in i|5,3l Although the groups overlap 
in their linear separation sensitivities, we quote a "Net" timescale 
that gives a limit on the total timescale of inspiral from galaxy 
virialization to SMBH separations beyond the stalling radius. 

is discussed below, followed by a brief consideration of the 
consequences of our results for the detection of gravitational 
waves from SMBH binary systems by pulsar timing arrays. 

7.1 Post-merger efficiency, radio ignition, and 
SMBH binary stalling 

The timescales measured for our three merger stages rep- 
resent the time a pair's orbital separation takes to evolve 
through the spatial sensitivity window of each stage. For 
our merger stage groups one and two (~virialisation down to 
SMBH binary formation) , it is expected that purely dynam- 
ical friction will act on the black holes. Formally, the dynam- 
ical friction timescale in equation [2] over the range probed 
by our first and second merger stage groups is ~3.7Gyr and 
900 Myr, respectively, for an object of mass 10 s Mq. 

For merger stage group one, the optimistic scenario 
timescale limit is smaller than the dynamical friction 
timescale by a factor of ~30. From this disagreement, it 
follows that one of the following must be true: 1) the rel- 
ative radio state of two SMBHs in merging galaxies is 
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not strongly correlated before the hosts' stellar cores have 
merged, and/or 2) The centralisation of supermassive ob- 
jects in a post-merger system is highly efficient (in excess of 
dynamical friction). 

Considering merger stage group two, in which we have 
detected one binary SMBH, the timescale under the opti- 
mistic set of assumptions gives a value more consistent (a 
factor of ~2 within the error bounds) with dynamical fric- 
tion. Because in the pessimistic limit we would not expect to 
see any sources, it is implicit that the pessimistic assump- 
tions are not strictly valid for this stage of merger. This 
is furthermore indicative of some radial dependence of the 
relative radio state of two SMBHs in a post-merger environ- 
ment, such that the SMBHs have a correlated radio active 
state at later stages of SMBH inspiral, after the pair has 
shared a common galactic environment for many orbits. 

We have set the cutoff for the third merger sta ge group 
to be sensitive to the theoretical stalling radius of Merritd 
(2006). While the limits that we can place on this stage of 
merger are the least stringent due to the small source count 
in the limiting group, they represent the first observational 
exploration of stalling. If the evidence exhibited by merger 
stage group 2 holds — that the relative radio state between 
two black holes is more strongly correlated at late stages 
of binary evolution — our results are in agreement with a 
progression through the stalling radius in much less than 
a Hubble time; that is, we find no evidence for systematic 
stalling of supermassive binary black hole systems. These 
results, however, exist with the significant caveat: if the most 
pessimistic assumptions are true, we would not expect to see 
any paired SMBHs regardless of the length of post-merger 
inspiral. 

As previously noted and reflected in table 2, we have 
only used the information from roughly half of our sources 
to acquire timescale measurements. The lack of detections 
in the additional 1539 sources suggests that the measured 
limits will only be more stringent, however would be highly 
dependent on the redshift distribution of the remaining 
sources. 



7.2 Consequences: gravitational wave background 
and pulsar timing arrays 

The measured inspiral timescales and lack of evidence for 
stalled supermassive systems are pertinent results for pul- 
sar timing arrays, which are sensitive to the most mas- 
sive of binary mergers such as those explored by this 
study. While we have primarily accentuated the brevity of 
timescale to pass through the stalling stage, we note also 
that the timescales measured for merger group two (in con- 
cert with time spent at other separations) are sufficiently 
large that they may have non-negligible effects on predic- 
tions of the gravitational wave background signature in the 
local Universe from SMBH binaries. Typically, estimates of 
expected binary black hole contributions assume efficient 
or in stantaneous inspiral into the gravit ational wave regime 
fe.g. ISesana. Vecchio. fc Colacind 12008 ), using merger rates 
as binary coalescence rates. The inspiral times measured 
here may cause the galaxies merging at a redshift of 0.7, 
at which the merger distribution peaks, to emit gravita- 
tional radiation at a much later epoch, causing a larger ex- 
pected amplitude for the astrophysical gravitational wave 



background than previously estimated. In addition, it is 
likely that this effect would cause an increase in the pre- 
dicted number of supermassive binary sources whose grav- 
itational radiation is strong enough to be resolved from 
the gravitational wave background by pulsar timing arrays 
(Scsa na. Vecchio. fc Volonteri|[2009l ). 



8 SUMMARY 

We have performed a search of 3114 active galactic nuclei 
for the presence of double SMBH systems using a multi- 
frequency radio imaging technique. Of the sources searched, 
only 0402+379 was apparent with our method as a binary 
AGN. While this s ource has already been put forward as a 
binary black hole in iRodriguez et all (|2006l ) , this search rep- 
resents a significant enough statistical sample with which to 
interpret the existence of this source in the broader cosmo- 
logical context of binary supermassive black holes. 

This searched has probed estimates of the inspiral 
timescale for binary black holes of mi , rri2 ^ 10 s Mq . For 
the 1575 sources that had redshift information with which 
to discern the linear scales probed in the host galaxy, we ob- 
tained timescale estimates for the inspiral of binary systems 
in three ranges of binary separation. We have demonstrated 
observational evidence against stalled binary black hole sys- 
tems by demonstrating that: 1) at late stages of SMBH pair 
inspiral, the two black holes are more likely to be in a simi- 
larly radio-active state, and 2) SMBH binaries proceed from 
separations of a ~ 5 0pars ecs to within the stalling radius 
estimates of iMerrittl i|2006T ) in less than ~0.5Gyr. The im- 
plications of this are that supermassive binary systems do 
not stall indefinitely at such radii, suggesting there is a yet 
undetermined mechanism by which the black holes are able 
to dispense energy to the surrounding environment and pro- 
ceed to coalescence. 

The results of this study come with the caveat that they 
stem on very small-number statistics; that is, we detected 
only one binary AGN in the analysis, and due to unmeasured 
redshifts, we were unable to use roughly half of the sources 
searched in the statistical analysis. Definitive studies of this 
type for large numbers of AGN will become possible with the 
large collecting area, dense instantaneous u — v coverage, and 
long baselines of the planned Square Kilometer Array. 
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